{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scalable GP Regression (w/ KISS-GP)\n",
    "\n",
    "## Introduction\n",
    "\n",
    "For 2-4D functions, SKI (or KISS-GP) can work very well out-of-the-box on larger datasets (100,000+ data points).\n",
    "Kernel interpolation for scalable structured Gaussian processes (KISS-GP) was introduced in this paper:\n",
    "http://proceedings.mlr.press/v37/wilson15.pdf\n",
    "\n",
    "One thing to watch out for with multidimensional SKI - you can't use as fine-grain of a grid. If you have a high dimensional problem, you may want to try one of the other scalable regression methods.\n",
    "\n",
    "This is the same as [the standard KISSGP 1D notebook](../04_Scalable_GP_Regression_1D/KISSGP_Regression_1D.ipynb), but applied to more dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import torch\n",
    "import gpytorch\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up train data\n",
    "\n",
    "Here we're learning a simple sin function - but in 2 dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We make an nxn grid of training points spaced every 1/(n-1) on [0,1]x[0,1]\n",
    "n = 40\n",
    "train_x = torch.zeros(pow(n, 2), 2)\n",
    "for i in range(n):\n",
    "    for j in range(n):\n",
    "        train_x[i * n + j][0] = float(i) / (n-1)\n",
    "        train_x[i * n + j][1] = float(j) / (n-1)\n",
    "# True function is sin( 2*pi*(x0+x1))\n",
    "train_y = torch.sin((train_x[:, 0] + train_x[:, 1]) * (2 * math.pi)) + torch.randn_like(train_x[:, 0]).mul(0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The model\n",
    "\n",
    "As with the 1D case, applying SKI to a multidimensional kernel is as simple as wrapping that kernel with a `GridInterpolationKernel`. You'll want to be sure to set `num_dims` though!\n",
    "\n",
    "SKI has only one hyperparameter that you need to worry about: the grid size. For 1D functions, a good starting place is to use as many grid points as training points. (Don't worry - the grid points are really cheap to use!). You can use the `gpytorch.utils.grid.choose_grid_size` helper to get a good starting point.\n",
    "\n",
    "If you want, you can also explicitly determine the grid bounds of the SKI approximation using the `grid_bounds` argument. However, it's easier if you don't use this argument - then GPyTorch automatically chooses the best bounds for you."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class GPRegressionModel(gpytorch.models.ExactGP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super(GPRegressionModel, self).__init__(train_x, train_y, likelihood)\n",
    "        \n",
    "        # SKI requires a grid size hyperparameter. This util can help with that\n",
    "        grid_size = gpytorch.utils.grid.choose_grid_size(train_x)\n",
    "        \n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = gpytorch.kernels.GridInterpolationKernel(\n",
    "            gpytorch.kernels.ScaleKernel(\n",
    "                gpytorch.kernels.RBFKernel(ard_num_dims=2),\n",
    "            ), grid_size=grid_size, num_dims=2\n",
    "        )\n",
    "\n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)\n",
    "\n",
    "    \n",
    "likelihood = gpytorch.likelihoods.GaussianLikelihood()\n",
    "model = GPRegressionModel(train_x, train_y, likelihood)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train the model hyperparameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 1/30 - Loss: 1.142\n",
      "Iter 2/30 - Loss: 1.087\n",
      "Iter 3/30 - Loss: 1.025\n",
      "Iter 4/30 - Loss: 0.960\n",
      "Iter 5/30 - Loss: 0.901\n",
      "Iter 6/30 - Loss: 0.843\n",
      "Iter 7/30 - Loss: 0.774\n",
      "Iter 8/30 - Loss: 0.686\n",
      "Iter 9/30 - Loss: 0.595\n",
      "Iter 10/30 - Loss: 0.521\n",
      "Iter 11/30 - Loss: 0.465\n",
      "Iter 12/30 - Loss: 0.412\n",
      "Iter 13/30 - Loss: 0.365\n",
      "Iter 14/30 - Loss: 0.317\n",
      "Iter 15/30 - Loss: 0.278\n",
      "Iter 16/30 - Loss: 0.229\n",
      "Iter 17/30 - Loss: 0.182\n",
      "Iter 18/30 - Loss: 0.148\n",
      "Iter 19/30 - Loss: 0.103\n",
      "Iter 20/30 - Loss: 0.057\n",
      "Iter 21/30 - Loss: -0.000\n",
      "Iter 22/30 - Loss: -0.060\n",
      "Iter 23/30 - Loss: -0.117\n",
      "Iter 24/30 - Loss: -0.178\n",
      "Iter 25/30 - Loss: -0.235\n",
      "Iter 26/30 - Loss: -0.290\n",
      "Iter 27/30 - Loss: -0.356\n",
      "Iter 28/30 - Loss: -0.405\n",
      "Iter 29/30 - Loss: -0.455\n",
      "Iter 30/30 - Loss: -0.508\n",
      "CPU times: user 30.4 s, sys: 237 ms, total: 30.7 s\n",
      "Wall time: 12.1 s\n"
     ]
    }
   ],
   "source": [
    "# Find optimal model hyperparameters\n",
    "model.train()\n",
    "likelihood.train()\n",
    "\n",
    "# Use the adam optimizer\n",
    "optimizer = torch.optim.Adam([\n",
    "    {'params': model.parameters()},  # Includes GaussianLikelihood parameters\n",
    "], lr=0.1)\n",
    "\n",
    "# \"Loss\" for GPs - the marginal log likelihood\n",
    "mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
    "\n",
    "def train():\n",
    "    training_iterations = 30\n",
    "    for i in range(training_iterations):\n",
    "        optimizer.zero_grad()\n",
    "        output = model(train_x)\n",
    "        loss = -mll(output, train_y)\n",
    "        loss.backward()\n",
    "        print('Iter %d/%d - Loss: %.3f' % (i + 1, training_iterations, loss.item()))\n",
    "        optimizer.step()\n",
    "\n",
    "%time train()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make predictions with the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOoAAADOCAYAAAAwqZ5uAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAFD5JREFUeJztnV1sHNd1x/9DiqTED2m5+rAoOZa88qdE1TY1SosEQb8oJH1oCxSShaIvfbH0VKR9iOgUaIOiaC35oUAf+iAG6HNtqwjSvqQxgwJum7TViHAayfKX2MipLcmSSEqk+KHl7u3D3pFG673nzu7dXc7s/H/AgMs5d+eenZ2z536e4ymlQAhJNl3rrQAhxA4NlZAUQEMlJAXQUAlJATRUQlIADZWQFJAoQ/U8b8zzvAue553yPO9o+LeB65zxPO+Ufl3wPO8toWys60evWXV+3PO8Oc/zTkTOnfI876zneTlD+bfjfZL4hJ8jcg/PeJ5XiMgf3Aetg3RPHuhY431N0b3qO8o18j1niUQZqlJqGsAMgCml1Dml1OsAvlvrgbfwRuSaM0qpY7UK6eseqfeaVTpPAZisOj2tlDqplJo3lP/CeRc8zxsHMK2vH97DN5RSM5F6H9wHrYORqI413tcs3aPf0bz+HAVz8WyTKEOVCH/NI542p/+Oh94s/B/AeOR9Y1EvoMuM6V/wAgBfvwdxr1mDswBORv5/8MMSXiuso9ZnitRzppYeWt/x8KhR/5GoURru3yP3IXL+gud5f1j9uaX36XKnPM8bi5wL9TxhOVfzfiqlzuHRe0giJNVQ/chD84pSal7/mhcinvbbqHjeKQCHdNlp/f8Dj6E9TPiLfRTAjD53PPQ+EQ8T65rVhEaim4k5VDxa6CEKSqlJABM13hf1UOciokf0AHA8Ur6WQVpbHNH7EKIN5hiA56vqk95X0OXOhXrpJmwQ3k9toLXO2e4nPaqBpBpqoJSaUkpN6l/akOnI6zEAef2rfhaVJqzoVaJlDM3hRq4ZEnrVcf1whs3GyQaa7tV6vAbgiOd5VxDDKKsxeOGC1jdXoz6JWvfjCB4a84z+33ROup+zlrozS1INNQ5vA4/0yc6j8sABQN7wniuhrNp49EPayDWh3zMJ4JEBEd3s+8IAlIGoN6nWY1wpNYGKt5Oa4CbGapwLALwC4EyN+uplGg/1L6By32qdi30/yaNsWG8FomhjKQA47nneTHQwRnuFMc/zxpRS00qp13V/BwAQ+X8MlYfhiOd5k/p6Y57nFXSZM/o9OVSabzO6STwV95q1Bok01U26AipepRCpZyaizwyA8/qz5QCMe56Xq9YDwOHI62gLIyR6n0Jdj+um9xEA8+H5sDkOwEfFWOYBbNWy8BqzkbI5w+vw+8gppSYi+o7prgmEc6b72dRBtk7C4+6Z9KMNfcY2oJRkOuEztJI0N32JRg/O1GrepoKwG0IjNUOPSkgKoEclJAXQUAlJATRUQlJA06ZnXn31VXZ2SWI4ffq0Zy8FnP6rP1LzCwO2YldPnz6911kpB5o6j/pWfodRtn+gH+/dWzLKj3/j353r/7Pt06L887kD2DF8yShfUWvOOnznxtdE+balXbjV/5lR/s/v+E71P/ZfsvypXf34+DPz9wAAm24WnXS4N9Ijyp8e6cdH12QdbjqMYR+/fjN22fmFAfzFN98Uy3znb1/e07g2zSFRCx4IWQ/KSH5jkIZKMk9RldZbBSs0VJJ56FEJSQFFlNdbBStWQ/V9/ygqi6XHgiB4vfUqEdJeSilYnSfOo2ojRRAEUwDmfd9vZIsVIYmmCCUeScDmUQ/jYWybGVQWfhsjHewf6Dde6ImNfWJFuaXHLarY+XxuVZTfWXxClBfhPqiwbWmXKN+8Km/DPLDJfA/jkJOrx0he/h4AoHfA7T6s5rtlHYbtOuzY5KRCXRSTYYsiNkOtjiawVSoszZPa5Af7/8+iih1pjjROmWbMo95aHbaXEeZRLy1bLM3CY+ZLP6Dl86ir8jwqAPs86kjj9Y/WWb6EWGsj1hWboc6DO/FJh1NU6TfU83joVQvQITsI6STS4FHFwaQgCM4BKISDSHpQiZCOoqi6xCMJWKdnOCVDOp1SCjaRccEDyTxJ8ZoSTTXU3/v6T4yy/NJuPNf/qVH+x1v/07n+FcsNL6Ikjuz+3eyLzjp8/39eEOWjGwZxcW27UZ674vbQqC55lY3yANtzeeugffpE4u4v3Rfl23vKuPG8PML+WwcvNq7AW7vrKl7KmqESkkaKkOd9kwANlWSeoqKhEpJ4OJhESAooqsbNQE9dTgRBUDN9Z7M2tST/p4SQFlNSnnhISGsLmrmphYZKMk9RbRAPBw7jYdKtcFNLQ7DpSzJPC/uodW1qkaChkszTwlHfpm1qoaGSzNNCQ23aphb2UUnmKaku8ZDQA0Z+OHCkz70NNHdTCz0qyTwuHlUb47mqc0cir5uyqYWGSjIPFzwQkgK4hJCQFFDm7hlCkg89KiEpIHOG+qfbzRvHb8+PYmvOvBl4tQmxVW0bv4fuPYGFsvkj//3PvuKsQ/+H8qbrvlwP+ueFMo7ZFeaelZtxi0NdmNssl+l9Yc5Jh28+bX4OAKB/cS9+c/DnYpnfGHi/4fr/AX9QV3luHCckBWTOoxKSRsqWHTJJ8Lc0VJJ5bB7VLYJUc6Chksxj86hJgIZKMo/No7YxX5UR0VB9388BCHelHw6CYKL1KhHSXtbKyR9MsvWTXwaQ1wuP4fv+idarREh7KcMTjyQgetQgCCYj/xYAnJXK3543J7xbWNwjKlJsQtbnoXty/tNNK+bA1wBwoHvQWYeenDz08PiALN9gz1gosjoky7+0yT400qPcJnP7F/eK8r6VHdZrLK21r1dWTIFHjXU3fN8vAJgNgmBGKictaLDJVx0fDgBYUPanfGHgE6PsUsk9mXLfvD0J8OV5c27Q3jtu9S/HeOYuL8i5SXs9NyW+YlnMAABLljL9DgseADlbQTWdNJh0NAiCky3VhJB1Yi0FK5OsGvq+fzTc/OoS7pCQpFJWXeKRBGyjvuMAzvi+/219iqO+pONIg0e1DSZNAdjXJl0IWRfWyik3VEKyQCcNJhHSsaS+6Vsvd8rmqYl75TI2CPK/uflrzvV//6fyftTRnkFcLJrncwc+6HXWYcOyLO/eCPTcM8sXnnSbT97zkjlZNABsXctj94ZZscy39v6Lkw6/3CfvZ71d7sbWze+JZabvWyaEmwg9KiEpgH1UQlIAPSohKcAlFIst/6nv+3OoZHKbctnUknyfT0iLKZW7xMNEzPynx4IgOOS684yGSjJPWXniIRAn/2lOr5V3goZKMk+jHhXx8p/mAcz6vi/uPLNBQyWZRyn5ELDmPw2CYDIIgnlUmsZHpbISHEwimcc+mGTcginmP9WBFmZ14IXbLjrSo5LM02gf1ZT/NMyPCuBNRAaZwkgpjUCPSjJPudz4PGqtKZkwP6pu8obJixtOYgzQUAmxDRglAhoqyTxNCNfVcmioJPOU6VEJST4pcKg0VEKUw2BSu6ChksyjsrZ75tQvfsco2726A58uPGmUnz//jHP9wx/IN3xgezeGb5pj/3atuTeC7lpWda4OKiz1muvZN/YLp/r/svA9uf47B9C35ZJY5st9blHA/2Nloyi/v9aDG0W5zJ9/9LsN1/+rdZZ3mZ5pF/SoJPOw6UtIGkjBaBINlWQeelRCUkDmBpMISSUpMNTYSzJ83z/TSkUIWTeU5UgAcdMujqOy346QzqMT+qg63ouYFzVk96o5Qe3W4hbxvUub+uNUIdK/Xb7huzfLSXy77KlNrSxZciF/aaOsw2PFWtE84rN654AoLy7JyZ4B4HqPW4/o/n051+1aDB2eKg07aDBfV+kmpOZtOXG+kUIQBFO+71sLftr3ecPyS8vV4WfqZ8tN+y/jhzfNSXy71pxVwN0YScvfWzTrsNzjFAjAupghTpmdjgserqzYn/zeLXKk/I9v7G+4/t2o00OmoI9qTbsY7lonpFPxOsCjzur+aQ6VkBNjQRBMt0EvQtpHCjyqOOobBMG09qh5fDE0IiGdQdlyJIBYowZBEEwCmGyxLoSsDwmZgpHgggeSebxOmJ4hpOPJmkd9d+o5o2x1qB+XF8xBxR+77N4Z6CrK1+jrUej/3Fxm7pluZx22HLwlygfKJWzpMs/z/cmet42yOLzUKy82+3xDF3ZYyryz4qQCTr3/sih/upTDR9dGxTIrPzDPyVvZeKOu4l7WDJWQVMKmLyEpwMGjxsiPKsrjkvw4iYS0GK8sHyZs+VFj5k+NBQ2VkMZ3z9jyo8bJnxoLNn1J5nGYnrHlR42TPzUWNFSSeRzW+tryo1rzp8aFhkpI44NJYn7UGPLYsI9KMk+jg0m2/KgmeSPQoxLiMD0j5Uc1yRuBhkoyD1cmEZIGaKiEJJ9OiPBASMdDQyUkDbDpS0jyyZxHHfnxfaMsv6cXI1fN8q6i+8/a3LO9onw152GpaJ46Xh1ddtbhlSf/W5T3L+7F4cGfG+WFnlmn+v9x8XFR3r2yGaXFbWKZv778DScdun8ox+QtbuvH6i35u9r1Izn0rEi9qVXpUQlJPpnzqISkEptHTcC+choqyTxWj+oeoccZGirJPFyZREgaYB+VkOTTER7V9/0x6NyoetsOIR1FGgw1zn7Uk5F9dUxmTDqPtOee8X3/BIALvu8X4uyrK+wxJwfduWOj+F5vzf1nbft2uYFgS2S83O3eE+hf3CvK+1bkwNKLa257+buL8mKGruUR6zWeU3LSaRtd2+Sk1LbvAQCGn3fJSTZXV+k0eFTbk7lP/33T9/2zACaCIDCGeZ+5uiheTJI3ZWVSn7zaBZATGd8tOYaIB/DrwqqjkCWhzOCgPRGxRGlFXpkEAKWhj0T5+94+UW6j+5b9x+aDW+bvAQBGLteXNTzK8DP1lU/Dgoc4P99XtHFeAHCixfoQ0n4aDxfaNmwe9TweRlHLoRJVjZCOIvUeVQ8i5SLBmZgjlXQcjQY3ayfW0ZPIIFLDEdQISTQJad5KcMEDyTxeOfmWSkMlmacTpmfqou+2eXqjJ98ryuf2DznXP/fimii/11vC3G5zmd9+9qKzDs/1XRPl5dUt6BLKfO/ui071f/dnXxXlB7qGcKm8XywzPLXJSYft/3Zdvv7zwxi5bJnrvCEnhLZoUFfppPRDJehRCcmaRyUkjbTSo8ZIdDyHSkrGqSAIJkzXYe4Zknm8shKPRomZyPhYEASHJCMFaKiEwFPy4UCcRMa5OJtdaKgk83gl+XAgTiLjPIBZvZbeCPuohDh4Tb3DrJqZsLkLSyLjcLWf7/vzvu8fNe35pqGSzOPSD7UsqxUTGWsjn9XGeVuqh01fknla1Ue1JToG8CYig0xSBBV6VJJ5Wjk9IyU61ttHwzX04lp6GirJPFzrS0gaSL6d0lAJ8UrJt1QaKiHJt1MaKiHsoxKSAjK3H3WhYI7ru/zYRiyUzOPgNw+7j5F/9aAcBnPX6g5s6TMnyPUH/9dZh39deF6U55d3Y9Yz5/F7891DTvXvfLtHlG8e6cH2a3KM5fw7nzjpUL5zV5SrnV1Q12+KZUoHXWK9y2Frq6FHJSQFcDCJkDSQfDuloRLCpi8haUDRUAlJPB3RR43EfCkwUj7pSJJvp/I2N739JtwEO6OTGhPSUXjlsngkAZtHDVDJj3oMFY8qbsV5apc5L+bIVjknZl6e2ovFrlU592i+KOfc7Fl4ylmH/P1a0TYeMrgqy0c3mOei45Ab6RblI3l7btKhUVlHG+V7clzgHQV77tPyXpc4z3XOo6a96RsEwbyO5fIWAOOm1pCPP5NzXkrya/vcb5a0mCHkM6HM6NDHzjrMLttzrM72f2qUXVzb6VT/zmvyggcA+Pia/D3lL4rBBqzYFjwAwNV3b4jyUmmg4fq32XM1P0oKBpNsTd+jqMQb3Rf5n5DOoqTkIwHYQrEUgiCY1q9fgyVQEyFpxFNKPJKArY86qQMwzYCjvqRTSciAkYS1jwqAxkk6m4R4TQkueCCZJ/WjvoRkAmH7ZVKgoRJibfqa9w+3i6Ya6vVfMcu2bgKuP2GW+y+5z2GO598T5T0Lq9g/dMUof395l7MOto3foxsGxbnSkR/Y50Elcj/8QJQPvrADwz+V55vX5ixJhi14hw7IBXZvhlqSF1588vXG51H9e/UteLAPJsmLSNoBPSohLdzmppfhToRBt2vIxfypIUxpQUi5JB8OSMtuY+ZPBUBDJaTiUaWjdcTJnwqATV9C1nPBQ5z8qQBoqIQApcabt5b8qDas+VNDaKiEOKxMclxWK+ZPjcI+KiGlsnw4oAeM/OjOszA/qil/ai3oUUnmUap1fVRtjOeqzh2JvDZOyUShoRLCJYSEpIC0b3MjJAsoh1HfdkFDJYRNX0JSQAsHk5oFDZVkHjZ9CUkBKmtJoo5bktOOSsJ/sgdltvEhvhajlPueU4nfx7ylxDwOSuI9jgq8YluRtobhL9vKuAabtMU2XkFe2JsMAGP37PGRm8TVQ9960nbXr7ZFEwFPpSCwEyFZh0sICUkBNFRCUgANlZAUQEMlJAXQUAlJATRUQlJAWxY8xA2J2KK6cwDC6G6HgyCYaGf9NfQ5s1466IzxBeDBPsl21x8+B0w4Vict96j1hERsES8DyIcPpiHGTVvQn72wXvUDOBmJKtBWPfRnD2MJzegfDRKTdnjUwwDe0K/DkIhxAj81hapf7gKAs+2qO4o2jBlrwdbVfwLABd/3C+1u1WgCXf8xVDxq256BTqAdfdTYIRFbiTaU2SAI1stYCutYNwDs08es7/tndZegbegUnmcBvKX1IHXQDkONHRKxxRwNguDkelTs+/54QjzIFW0wFwC0tQugu0BTQRDsi/xPYtIOQ40dErFV+L5/NGzurUMfGah4sXH9cBbWqX92PvI6B1h3DzSbQhAE0/r1a0jGj3dqaMuifN/3TwGYxjqM9mnDPIuHD+bEenk33U+cAHAs8tC2s/7we1iv0feXUemnc9S3Trh7hpAUwAUPhKQAGiohKYCGSkgKoKESkgJoqISkABoqISmAhkpICvh/KqONs950soYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAO8AAADOCAYAAADWgFUqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAEixJREFUeJztnc1vJUe5xp+aQUkmH+SMIUMAQcIxLPhQLnjKGzawsPkL7JkViyywFxG6EhJ2wgbdxVVmsrvL8ZLdzHjD9uL/YGq8uAiQUMYiC5QoEM+BTCCRx6fuoqt9etrd9dbpOl/V/fykVs7pqu4qT/z6ra/3eZW1FoSQ9Lgw7w4QQppB4yUkUWi8hCQKjZeQRKHxEpIoNF5CEiU541VK3VVK9YQ6Gw3fvaaU+l3N/YdKqa3CvR2l1K2qvtS9J5b851JKrSil7iulbiql+oXyvlLqbqEPd4X3nfWz4tmJ9N/1ccd97jX9f0POk5zxAugD2KordMa03uTF1toDAIOa+3ul24fW2m1rbV39c/djUEqtATh07z8EcATgtrX2qNDukbV2s9AHL8V+Vjw7qf7fLrQ3cD9Lv746CSUp43V/tX8G4Hrp/o7zRhvIjFu7X/ayd9lRSt0sPLemlNrK6wrcArBd+H7mcX3vqWvfeaGdwrMr7vNaTX/Wi4ZahXtH1cjhvlJqo9ym9Kyrt6OUWincy/u6JdzbcT/HEz+LtXYfT/47koYkZbwAlpzX6eV/vZ3BHrn713OvlHuekhfZz1/knu9ba/cA7EoN54bjhpc9ZJ5PfE9d+wDeAnDgyq/C/UFy36uM1DtVcM8eouQxnQFtOqMpt+l7tu/q7ed9c8Nfk/8bO6OtureFbGRyAKBqBEDPOwFSM96r7pfxCEA+d1p335EP+0Jww8Q9af5cIve+a+6Xtel7AGAFwJLzarcAvA1gXSn1AAGGWqbGW/ddf/P3ldv0UfUHZB0jAz9y3+vu+UYJx0LbJIBkjFcpteHmmAfIfiHzofMDAEuuTq/0zMqTb0FxcWclX0gJxXnXJxZcxnxP0eP8zr0zn7+uWWt3kXnEkGF8mfLPCgAG2TQjnyqU2xyXQ4x+hj6AezX37hX6s9SgHRLA5+bdgRCcEW4rpQ7d8LUHoK+U2rLWvuNWNOHu7yMbvm1gNGS75zxTD8CaM/I+Mo/RL9Q/ArCilOp75pfloWDIe8617/q94/oNAKuFz8Xhdc7ZkNb9e6wAuO6G7esABvn9fCgPQCMznoFS6qa1drfYplLquFC/V/N5zX3ulZ5fsda+495Tdy/v57pSaq+wuDfRxbyuohhVlAb5dEFatFp02vJzLALJDJu7jpsuVA2NkyGf1tBwJwM9LyGJQs9LSKLQeAlJFBovIYkysa2iN998k5NnsjDcuHFDybWAG//9czv4+Dmp2ns3btx4NbpTE2ai+7y3X36ptuy7l57FH/79r9ryN37yv9Ht/2LJv4j5wcPv4OXLf6wt//PJJ9F9eP1PP/WWf/P0Mt69+LC2/LP9L0W1f+W373rLX3ntJbz3f3/z1jn9m79cwv7wP7zlry6/gL88+Nhb58HmpcbtX38/vP+Dj5/Df/3nHW+dX//PtVcad2aKJHFIg5BpMkSag0YaL+k8J/Z03l1oBI2XdB56XkIS5QTDeXehEaLxaq03kB0kXzHGvDP9LhEyW04TPWXo3ed1hgtjzAGAgda6SagaIQvNCaz3WlQkz7uKkQbREbKD8bXaSN+99Gzti77+9NPehp57FL8a/4F6xls+ePQ1b/mjx59G9+Gbp5e95V8ePu8tP7lS/28Ywouv1W/XAcCVV18U3zH8x1NRfbDLL3jLv/SyvA30zKW4PozDyeLapxfJeMuKDl/wVfbt40rlP37+PaErMi9floNVfPu8/5zAPu+7739fruPb5/0w7pf2irCHC2D6+7xPfVGsI+7zrjTf5/3emPVPEXSeY+GQjHcAKiGQlnNi22m89zDyvn04GRVC2kSqnte7YGWM2QfQzxeq3MIVIa3ixF7wXouKuFXE7SHSdk4TDa7jIQ3SeRbZu/qYqPH6IoOee/SKd0VZiggKQYoKevT4U++KshQRFIIUFXRy5VnvirIUFSQhrRQP//GUWEeKCpKQIoKeufSUuJr8xnrzKLO//+YHY9U/pfESkiYnuDjvLjSCxks6z4ltbry+48Na6xUAdzHSqT4wxuxqrR8iO/R0YIwRU+3UQeMlnafpglXx+LDWuq+1XivtyCwZY5Zd3RWMjHhzEjs3aQ72CZkgJ/Zz3svDKkZpY/Ljw2eUDLRvjMnr9rTW0cnWaLyk85xa5b08BB0f1lpvuTMTOUsAjrXWUrI3LzRe0nkiPG/o8eEnkr0bY/aMMQNkkXobNc+IcM5LOk/EIQ3x+LDWulf6vgXg2Hnij5o2DNDzEoITe9F71VF3fFhrXTTiJTyZj/gOCrHxpeH0WNDzks4Ts1VUdXzYGLNe+HyELJ90/n2AUUx81IozjZd0Hp6wIiRRYjzvPKHxks7DqCJCEoWel5BEGXLOS0ia0PMSkig0XvgD6j9Qz3ilWWeWXtMjzRqbXhOQg+lffO0lrzzrtNNr2uUXRGnWmPSagBxILwkzAHHiDL8Cg/EJ6QT0vIQkylDQbV5Uv0zjJZ1H8rz+RD3zg8ZLOo/keRcVGi/pPJLnjVu+mx5e43WxiHlaz9UYsSxCFpXHw+kI0Lnyc2Jzk8p5Lc3FryET0dp3jW41bYiQRWUI5b3qCMxfvWmMuVoy3InkvPZ6XmPMXuFrH4BXc+eDh9+pLetCblwgPj9uG3LjSrmWn/70ivgOKdfyJDlp7nlD8lf3tNZF8bmxcl77CJrzOqW740IHKvHlvpXK25AbF4jPj5t6blwgLNfyJ0KdkFzL9Xx7rNoRC1YhAnRnYnPGmO3AZ4IIXbDacA0T0joeNz9hJQrQ5aNXrfWgMNedSM5rsdda6418Uh0zPidkURnaC97Lg1eATmu9VVCH/CjkmXHw9swZ602t9X2t9f2mjRCyyDy2F7xXHQECdOfE5iaZ81pasDoAsNz05YSkwONh8wOQPgG6OrG5SeW85iEN0nl4woqQRIlYsJorEzVeX0xuFxJbA/HJrVNPbA3IsbhSbDcwmfjuUOh5CUmUmDnvPKHxks5Dz0tIolAGh5BEOeWwmZA04bCZkESRPe9wJv0YFxov6TzWzrsHzaDxks4jL1jR8xKykHDOS0iiDIc0XkKSJGaryCcmVyfgWCVK14Q0N7gImSDW+q86AsTk6gQcnxClawqNl3Se4fCC9/KwisyDAiMxuTOMMXsFEcc+RjG9PacLFwWNl3QeK1wegsTkKgQcz0TpGncaNF5CYIfKe3kIFZN7QsDReeQBsqH2huc5L1ywIp3HNt8qEsXkKgQccy+8j5EoXSMmary+gPouJLYG4pNbp57YGpAD6SVhBiBOnOFHY9ZvulVkjNnXWu9UCdAZY9YLAo5vuUd2kYnS6aIoXaPGQc9LiDQ09iII0NUJOJ4TpWsCjZcQnm0mJE1iPO88ofGSzhOxYDVXaLyEJGq8wfu8Wuub0+wIIXMj4pTGPAlN8ZnvTxHSPto653VHu4I2/nzJrbuQ2BqIT27dhsTWUqL0f3/i7yMgJ0r3Mxirtl3MWHuREM/bN8YcaK3Fir7E1VJ5GxJbA/HJrdNPbB2WKP353u+95SGJ0uv4Ksb0pInOeb3Gq7Vei0lBSEgKqJZ63mM33+0hyym6Yow5nEG/CJkdiXpe72qzMebQed4lnA9/IqQdDIVrQQlabXYBxXtiRUJSZIG3g3zwkAbpPCpiq8inYVVXLj0TCoPxCWl4SEPSsKoqD9C9CmaintcXk9uFxNZAfHLrNiS2lmJxpdhuIDK+e+nDsaqr5sPmVQC33edcw+pAKP+C8EwwHDYT0nzYLGlYVZUH6V6FQOMlpLnnlTSsqspDda9EaLyk80Qc0pA0rKrKe8IzwXDBipCGC1ZOf6pfpWFVV173TBPoeUnnidkq8mlYecobbw8VofGSztPWs82EtB+esCIkTeh5CUkVel5C0iTihNVcofESQuMlJE045yUkUWi8hKQKh82EpAk9L/wxuV3IjQvE58dNPTcuIMfiSrHdQGR89+ufH68+PS8haULPS0iqSJ53QZVhabyk84ie9+JMujE2NF7SeaZ5wsqnFKm17gHIBehWjTG77v5DZPpWB/m9KhiMT8iURNcDlCKvAVhyAfrQWm+5+5vGmKs+wwVovIRAWf8VwSpGGTZzpcgzjDF7LqEBkEni5KoaPZed04tovFrrFa31Rv5XhJC2MUXjDVKKdIZ6bIzJDX0JWZ6wW76Xh8x5t40x21rrHa11v9AAIe0gbmi8VXH7KB8qI0wpcsMYs51/yb2x1nqgtd7Ih9VlpBSfWwDuO6MVdXde8SS37kJiayA+uXX6ia3lROlfefFp8R1SonQ/n41VO8a7Foa9VUjqknDGmadBWXP1jp3BfuRrW/K8y+6/d5wL3zXG1KYd9yWulsrbkNgaiE9unXpiayAsUfqfP/yXtzwkUXodl6+Od8JqWoc0jDH7bsR6Tl3SGLPu7t/UWr/lHtkFcCercvZMpdcFwobND4wxA631fQBbACaifEfIwjDFrSKfuqQz5uVzD40WrryysJLx3sNozN5DNoYnpFWkejxSSq69j2zZOnfhzNFLWoca+q9FRRw2F9x+Y2V3QhYaRhURkiZqmKb10nhJ56F6JPzbPV1IbA3EJ7dOPrE15EB6SZgBiN06XIytomlDz0sIPS8haULPS0iicMGKkEThghUhiaJO592DZtB4CaHnJSRNpjnn9WlYufJzelXSMzmUwSGdZ1pKGgEaVkBJryrwGQA0XkKmGZjg1bBylPWqQp4BQOMlBGpovVcEIRpWZb2qIN0rgHNeQqIWrGI1rMp6VSHP5NB4SedRp82tN0bDyhl+Wa9K1L3K4bCZECtcDXFG2a/SsHJV7qCwKGWM2a97pgp6XtJ5prlVJGhYDVChVxWi1ArQeAnh8UjAH5PbhcTWQHxy6+QTW0OOxZViu4HY+O5HY9VmYAIhiRKzYDVPaLyEpGm7NF5COGwmJFUsjZeQJGntnLdwZKvPjAmklaRpu/4TVu6UR35O80hrXRvhQEiqqOHQey0qkuc1yPLzbiLzvN6UJ6968uN2ITcuEJ8fN/3cuHKu5Svfqg2UOUPKtexnzH3eNg6bXWrPWwDuAqjNE5rjy30rlbchNy4Qnx839dy4QJhg+l/MX73lIbmW6/jil8d8INEFK2nYvIFMnmO58J2QdnFq/deCIkUV9Y0xh+7z2wiMMyQkJZS13mtRkea8ey7m8AhcbSZtZYqLUj4xObcAfBejpPUHxpjdKlG6KsQ5LwAaLGk3U/KuRTE5rXVfa71WWvRdKkxJVzAy4k1pcRhgMD4hUKfWe0XgFZMrGWjfGJPXLYvSVULjJeR06L+aEyQmp7XecgoaOWVRukp4PJIQcdisaktiBegc6yhMT8uidCXDPmOixusLqO9CYmsgPrl1+omt5UB6SZgBiBNn0O+Pd0hDXrC6WFsSI0AHAFrrXul7lShdJRw2EzK0/qshAQJ0gBsiF76fE6Wrez+HzYQMpTSBzc3EJ0DnPh8B2C58rxSlm2yvCGkLDMYnJFEWOHLIB42XkNM0s2vTeAlZ4PPLPmi8hMQdxJgbNF7Seayl8RKSJvS8hCQKV5sJSRPL1WZCEoXDZkIShQtWhKQJh82EJIrl2Wbg+vv+ONDvecr+/psfRLf/K4S849vR7fj4kVhjgK96grux9GFcB17/vFDhM1y+KtWRyiWkeNpHorby2DG5zXnv6i+/4VfCB2RB8DmgbKJHwwjpOgzGJyRRaLyEJAqNl5BEofESkig0XkIShcZLSKLM5JCGL9nSDNruAVhzX1d9iZtm1J+b8+qDy4fTB/ySolNsP/89YNK6CTB1z1tMtoSCHu0MuYYsodO+60+Vwv1McD+7mINmimwXtIRn2g/3s+eZBI7cHxISwSw87yqA2+5znmxJzIA2KUp/4fsAvPlfpoUzFn9Kh+m2vwXgvta6P+vRj8O49jeRed6Z/Q60lVnMeYOSLU0bZzzHhUxss6Y/x7YBYNldx1rrW+U0G9PGiYnfQpaPdnmWbbeVWRhvaLKlabNhjNmWq02eirys8+KBM6L7AGY6fXDTp4NCPtqNWbbfRmZhvGKypWnjMq294z7Pes4NZN5uzf3C9uc037tX+NzDKJHzrOgbYw7d57exGH/Qk2YmgQla6x0Ah5jDKqMz1lsY/bLuzssLunnnLrLM54dS/Sm0n/9/mNeq/zVk836uNk8ARhURkig8pEFIotB4CUkUGi8hiULjJSRRaLyEJAqNl5BEofESkij/D1R+F3h6TjJ9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOIAAADPCAYAAADoIg0/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAFARJREFUeJztnc9zG0d2x78D/v4hEiT1W5YsQ3FtnI3XCdRM5Q+AKqfcaOWwlWOoU7KnlZxTKoeUJd9ypP4DW/oPzD2lanPQUKmtrKK1d0lLsq3VyjJFU5YoET86B/RIIwj9eogGMQPM91M1RQDdmH4DzpvX3a/7vUBrDUJIuhTSFoAQQkUkJBNQEQnJAFREQjIAFZGQDEBFJCJBEJTSliEPDJwiBkFwLQiCYux9JQiCax2eqxIEwecJ6i0lPF85CIK1IAgumnMvBUGw1olsCduLt3Oxg++vACgmvT7SOQOniABKAJajN1rr1U5PZL67JdUxSn8u4fluAtgAsKq1XtVaXwdwKf7g6BZGeW7G2lno4DTzWuub5vtkHxkoRTQ33z8B+IeWolLMMlRM3bL5rBL7bNl8vtzy/deso7E0V6JzA1CxcxRjluiN88TOVwyCYMkoylZ0/sh6tb5vJ1+7OjE20FTyqGv5sXQdbdovR79brHw5eh/7ftl8J9F1Ewta64E5ACybv+sASrHPP4+9XjN/rwComNclABcBlM37Suxc12LfvRarf6X189h5o/OstJHxmmnrYtRGrGzd9l6Qb721jfjvAWDN/B7lBNfR2n68XtTe5+bvEoCl2DWJ181DPgbKIgI4a57YG2jeKBHx7uWmsRIfAzgXBME6gKh7GdXbQMLuZhvKAOaNRVmx1FnVWn8CYBVoWmfz+c2WevH3NvlavwNzzqLW+qrW+qype6VdPaG9l2itN7TWV1u60OeMHNBaf4hk100sDIwimm7eBd0c113A693T+A20pbXeQNMaXgJwFk0LcxPNJz/M3xtCc21nEs1N+Dnw2njQitZ6w9zc81I9w17kA4DzkeKY62031k00I2q6n61d33UYuU07ia+bvMlAKKJRgAux8VARzfFNNFbZCIKg1PK0XjRjygqA60Ypy8ailrXWn5j65dh5b5jyIoBKzEJsmHNtGEsXnUe1kbNsvlsx3/mV+X7FfK9s6r723iLfa3XaELWzjNet1BvX0ab9cux9CU1FLkXXaq7znPleRbpu4iYwfXpCSIoMhEUkpN+hIhKSAaiIhGQAKiIhGYCKSEgGGE5bAEL6GaXUEpqunXIYhp+0lJXRXHUU+XBXwzC81O48XVPEjz76iH4QkhkuX74cJKr3H/+st55MuardvXz58unWD40SIgzDVaVUSSlVCcMwvslgPgzDM6ZuGcIGgq5axLVP1q1lp9UJ3Am/tZY///u/8W7/YVm+nPdmJnF7+5ldhpO73jIMTdbE8p8WDuBW44m1vLE75CfADyNy+xOTuLVj/w0AYOKB34hl7su6WH7mxBTWv30q1pn99d2O2//gH52K9ZKtJ1P49198Jtb5t/88/7alaBHAp+b1BpqLNV4qYotSlsIwtO5iYdeU5J4GOu7MtW5fa7vVTCm1HIbhVelEnKwhuaeq6+IhsIVk64SdGwhoEUnu8bCIN/DKKpZgFr7HUUol2vTttIhKqSWlVEUptedQC4T0A1U0xMOGGfOVlFIV834VAJRScYWcB7DpkkG0iAlmhQjpe+oeGx9aXRbms3Ox1xtobssTcVnERbzaWxbNChEyUFShxaMXuMaIiWaFIk6rE9ayw+/KsYt2TyafcraxMCNP/Z+cHBPLqyOj3jIUCvLU/anChFjeGPacP5uQf4NTo/JvAABjRT8Zpk7Yu3MAcGzBLcPU+wc9JNjZU+1qBjzgLkVMOisEAKKf0FX+/JhdiZPy8JB77kn0I852wY84IvsRAch+xJqnH3FH9iMCcPsRtzz9iN/KDyMAbj/i/z7quP3iX+3toV5HIt//vuK6c52zQoT0O1WdviKKjz7brBAhg0QdgXj0Amdfrt2sECGDRFWnv66FDn2Se+oZWGBGRSS5Z+As4u7f2aPo1d6exu7CMWv5D+/4i/L8lDzrWR0ZxfOivc5fv9v5iv+IcvFrsXzm6SksTt2zlle136zpjU3bRoEmb1UXUB/5Xqxze+wtLxkmH8jXUJsM8GJGvvkbCz7pQKp7ql0fNEUkpB+pwtNl1AWoiCT3+PZCugEVkeQeTtYQkgGqOn01SF8CQlKmnoGVNVREkntoEQnJAD5jRCmcoimPsmlBCh6V/iiVkJSp6iHxsBHfOA9gK1qT3cKF2Jptaz5KWkSSezzcF2I4RaXUMoA1pVTJtWabFpHknrouiIeAa+P8GXNsKqVWpEBSVESSezrtmiLZxvn1MAy3AKwBWLZVYteU5B6PyRrXxvkbeKWoRQgh92kRSe7p1CK6wima8mKs3BrtmxaR5J6Gx+6LBOEUo3IxugUVkeQeLvomJAMMnCI+et8eF/Tg3DAezdjLn55yh+Bz4drYe3L3EEZGv7OW/8uJX/nLMCaHCfy+UMVC8bfW8vs1vyCb1YZ8U809A+YmvxHr/G7yqJcMtUn5tqqPBahNyus7n791wEMCZ4T71+XhxmBC0mfgLCIh/UjDsfuiF/aSikhyj8siuhME+ENFJLnHZRF7ARWR5B6XRZTTBnUHV37EIoBoa8diGIaX9l8kQnpLzTHT3Atc49DzAOajDY1mWwchA0UDgXj0AtEitqyNKwFYker/+dykteytaXnI+3xczqmXhJO7h8Tyg7VZsby2/Z63DN+PPBfLn/woBwB+6ulOnXsmBweefiHnqQSAnw5Ne8kwdVDOM3liJkF+xJo7vZ2dvfkRXb7XXpBojGh2Fm+aNMRWfvdYzrsnlT894O/Ql5z1EV8LdYZnbnvLsOBw6AMQHfovPB36jxvu3ICPHQ79W3U/h/7so3FnnS8eyfdK8e7eonXHOfKzvdXvp8mapTAMnXnACelHahlYWeOUQCm1FK0gt8TkIKSvaeiCePQC16xpBcAVpdS/mo84a0oGjixYRNdkzSqaMTcIGVhqjX0Np/gYzcBSq5L7L/1HASEp09CBeNhIGE7xwzAMz7p88FREkntquiAeAotoWjvgVTjFVopSPNOIri5xe7Fgn3qvTmu8KNjLhw6+8G7//dn7Yvns0yEUp+x1SiPb3jLMFmQf3E4whtmCfdHUfcjT+i6eNWQf3mRj2FlH1/2m8xuOu6ox5K6zO9M7356H+8IVThFoBo/aVEqtSJ4HWkSSe2qNgngIOMMphmF41YRT3Iq6su3gom+SezwsohhO0SwJ3TRLRMV86bSIJPd0GunbFU4RwGeITeJISWhoEUnuqXu4L6RwiqZLGoVRZDhFQiT6aa0pIQOL2yL67wxyQUUkuUf7bXjpClREknvccU1pEQnZdzhGJCQDNBpUREJSx8d90S2oiCT3cLKGkAzQoEUkJH0yYBCpiIRoTtYQkj560NwXhV37BQXVQCyvduGpdG9H3BqGo7vTeFCw1/mvcTn4bxKmCo4Nzs/Hgaf2zcO/efYXXu3/+sE7Yvm7jSJ+vy0H+B1+KG8cdjH6RO7sjUxqZ53Rbf84t0mh+4KQDMCuKSFZIAOzNenP2xKSMroRiIeEUmpJKVVRSl101LsilVMRSe7ROhAPGwnDKUaBusVIblREQnQgH3ac4RRNKEUxeROwB0V0mVZC+hbtOOwkCadYcmVRA5KnZXOaVkL6ls5nTcVwikqpShRQyoVTEZOaVgB474A9UenJCdl3VR/yDyh79PkRsXxut/UB9jrDT3a9ZUDgyOv37IRYPPtCTrbq4t2GfI3HE+RP/HHK/n9MwrT8b8DxojtR6fiIjx9xb4lKdef7fsVwimgGFq6YOiWlVDkMw5vtTpTEIpbCMFxVSjkr3n4iR6mWyquHdxKIIrMw/idnnQdCndqBP3jLMOZy6APAzJfWoh+eyRmHXfz+8Yy7TmFLLH/81M+rVXT/G/Dln+R7Zfp+5xmD95qo1DEOtBKG4XWl1MV24RTDMDwXKZ2Jbyo+IZ1p2ZKaVkL6lcAjEoYUTjH2/iqAq9J5XI++xKaVkL4lA2tNxVnTMAxvGos4D4dpJaRvaTiOHpBoMJDEtBLSt2RgiRvXmpLcE3DRNyEZYNAs4riQeGrUUb4757cHDgD++95psfy9YAa3tT1J6K3vj3rL8KIq/6Q/0bP4IviJtXx70+3nkxi/J/+OOzOTeLI9ItaZ3/C7Mw/ck104k3oEs1/JdcbWH3YuwM/cfso4waApIiF9CbumhGQAWkRC0sfHod8tqIiE0CISkj50XxCSAdg1JSQLsGtKSPrQIhKSBWgRCUkfn5U1JpLbFoByu72Jschu58IwvGQ7D6O4EdJh8ChXOEWlVBlNBVwFUDZhZ9pCi0hyj8cYcRHAp+Z1FE7xZUQLs4n+plKqCGBDiuZGRSS5x0MRk4RTBACFZvfVCrumhHQe11QMpxhhuqbFqCvbDioiyT1BQz4ExHCKSqkrJoIb4FDarnZNJ76zSz063BDLq1P+cU0b3x8Qy4OZCQTb9kfc0135+0kY3ZbLCwcnUXhkf/4df+iXF3Dqm6di+cEzBZxcl+uM3N9bXNBWGt8JG08BDFWPYuR/Hoh19AF7Dkk3e9uP2Kn7whVOEcAKmkHXKgCKJuRMWzhGJLlnv8IpmsmZaIJGDEtKRSTEZRF7sCacikhyj9Mi+o+anFARSe5hzBpCskA/LPo2y3RKQHOWaN8lIqTHZMEiJvEjXjAKWJLWyhHSrwRaPnqBKxvUMoA1pVSp3TQtIQNBH3RNz5i/nymlVgBcCsPQumbuz47ZE1wem5edrDtF/0U+DcfVnJyUZSh0npLvJcOOOMknZmQZxkb87oqxMTl48JFj9gDLEcNzfvP1jR9kGQ6X5pznCCbcctrZW47JLHRNk0zWrIdhuKWUWgOwDMBqGf/wRzn5pFS+PdYFRUwQLPz2tl2GQhcSBrtW1gDAF4/sMkx6r6xx34R31p+I5fu9sgYA7jpW1gQeK2vmyntbIZWFHfquu/9G7HURjhXkhPQlnS/67hqu/IjX0Vw1Hq2lY2o2MnB4LPruGs6uaWyShim8yUCSha4pHfqE9MlkDSEDTdBIXxOpiCT39Iv7IjFj2/ap95FiQyyfuePf/ugTubM/e7KOw1/bnYXjD3e8ZRh2TP3P/+UCjv/WPr1f++Zbr/aHFuTIDUMThzH8pSMJ6KTdH5yEwpFDcnlxFoUjspumdmTWQwJ543MrHCMSkgX2Ka6pid4WhVhcZFxTQgQ6dV+44poCOA9gPtosEYtf8wa0iCT3eEzWuOKaxv3uJTRj2LSFikhyj8dkTaK4pmbX0iYDDBMiEHS+vDdRXFMAS2EYXpAqcIxISOdrTcW4pkBzHBlN4rQZQ76EikhyT9DQ4mEjtmH+jbim5m8FwBWl1JrZvWSFXVOSe3wc+o64pqt4tadXhIpIcg8d+oRkAK41JSQLpK+HVERCgnr6mkhFJCR9PaQiEsIxIiEZYOD2I45uVe0NzdfF8olvfvRuP7gnh+gb/+Awpn5j34und/z3I+LYEbl8aAgYsf/sweL7Xs3vTssxRWunp7Gr/fYbuijsyv6A2vw0qsflu1/3IBVaBC0iIRmAkzWEZIH09ZCKSAi7poRkAU1FJCR1+mKMGAuOU2LIfTKQpK+H8n5Es59qw2zn2DDZgwkZKIJGQzx6gcsihmgmKv0QTYso5r84XbKnwzpyVM53V3g+7hDFTTAn/2iH35FjZepd/7xshbnWMCYtMpyS0401Zvx8fPWJIbH8yBGfvIPJCKqyiXHdC82T+Eggp517oymPrqkUTtGUV9DMK3pOOo+oiCYv4gqAawCuu4S6syH/AFL50I8vXKd3EtxzBM4FcHefHfqFY+4n6N3/swchrh/0y5Zaczj0AeCrr/wXT0i4HPqA+17xcegfPLrHL3Q4WRMPp6iUKimlKq3GypRZ45lGuLqmSwBWwzA8E2+YkIGiruXDziKaYRSBV+EUO8IVs6YUhuFN8/pjJItYRUhfEWgtHgKJwikmwTVGvGqiE2+As6ZkUOl8QiZpOEUnzjEiACofGWw6d+g7wykmheEUSe4J6lo8bLjCKZrXS80/8vwKV9YQUu/cVyiFUzSvryOBx4GKSIiza7r/myO7qogjd+w+uuHpOkbuCAk67//Ru/2hA/YFBQCag/K6PdFB4dBBbxnqC7IMjZkJsc6Op8N9Z0F26D87NILt+phYR3sOWMa3ZAuzWxzBzmFZhgO3HvkJsReckzXyb9oNaBEJ4TYoQjJAw5UOav/VhIpICC0iIRmgRzssJKiIhAgTeL2CikgIQ2UQkgE8HPrdgopIco/WVERC0ocWkZAMwFlTQtJHc9aUkAzArikhGYCTNYSkD7umhGQAPWhrTT/4ubSX7hmKH0jlpW6KYqGBub89tM9tuGKj7mD+pM/3PdHA245cqt44f+JHOP6uo4qrvHvcPfvLd9521dlvIQKdgeU9hOQdBo8iJANQEQnJAFREQjIAFZGQDEBFJCQD9MSP6Moht89tFwFUzNvFMAydKbL2WZ4raclgEs2WgJeBb3vdPrNPW9h3ixjPIQdgKwpP3kPOA5iPbjyTVCcVzLX3wmFq40IsTHxP5WD2aZleWMRFAJ+a11EOOTHzcDdpefKWAKz0qu045sbfcFbcv/aX0cz+XOp1r8Swp+zTeaMXY8Su5ZDzwSjCZhiGaSlDKcW2AeCMOTaVUiumy94zTGaxKPv0mV623Q/0QhG7lkPOk6UwDC+k0XC7lM4psW4UYg1AT7vozD4t0wtF7FoOuU5RSi1F3bEUxqhA0wpVzM1XSml8dCP2uojmA7KXMPu0QE/WmiqlLgK4iRRmy4zireDVjXcpLetkxmmXAHwYuyl72X70f0hr9vo8mH26LVz0TUgGoEOfkAxARSQkA1ARCckAVERCMgAVkZAMQEUkJANQEQnJAFREQjLA/wMcMV12qHuLpAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Set model and likelihood into evaluation mode\n",
    "model.eval()\n",
    "likelihood.eval()\n",
    "\n",
    "# Generate nxn grid of test points spaced on a grid of size 1/(n-1) in [0,1]x[0,1]\n",
    "n = 10\n",
    "test_x = torch.zeros(int(pow(n, 2)), 2)\n",
    "for i in range(n):\n",
    "    for j in range(n):\n",
    "        test_x[i * n + j][0] = float(i) / (n-1)\n",
    "        test_x[i * n + j][1] = float(j) / (n-1)\n",
    "\n",
    "with torch.no_grad(), gpytorch.settings.fast_pred_var():\n",
    "    observed_pred = likelihood(model(test_x))\n",
    "    pred_labels = observed_pred.mean.view(n, n)\n",
    "\n",
    "# Calc abosolute error\n",
    "test_y_actual = torch.sin(((test_x[:, 0] + test_x[:, 1]) * (2 * math.pi))).view(n, n)\n",
    "delta_y = torch.abs(pred_labels - test_y_actual).detach().numpy()\n",
    "\n",
    "# Define a plotting function\n",
    "def ax_plot(f, ax, y_labels, title):\n",
    "    im = ax.imshow(y_labels)\n",
    "    ax.set_title(title)\n",
    "    f.colorbar(im)\n",
    "\n",
    "# Plot our predictive means\n",
    "f, observed_ax = plt.subplots(1, 1, figsize=(4, 3))\n",
    "ax_plot(f, observed_ax, pred_labels, 'Predicted Values (Likelihood)')\n",
    "\n",
    "# Plot the true values\n",
    "f, observed_ax2 = plt.subplots(1, 1, figsize=(4, 3))\n",
    "ax_plot(f, observed_ax2, test_y_actual, 'Actual Values (Likelihood)')\n",
    "\n",
    "# Plot the absolute errors\n",
    "f, observed_ax3 = plt.subplots(1, 1, figsize=(4, 3))\n",
    "ax_plot(f, observed_ax3, delta_y, 'Absolute Error Surface')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
